﻿#ifndef Shell_H_
#define Shell_H_
#include "StructuralElement.h"
#include "../../mat2d.h"
/**
 * \brief 壳单元的基类
 * 壳单元的刚度矩阵分为三部分平面应力、弯曲部分和剪切部分
 * 参考文献：A Smoothed Finite ElementMethod (SFEM) for Linear and Geometrically Nonlinear Analysis of Plates and Shells
 * 默认是六自由度
 * 应变张量为{sigma11, sigma22, sigma12, ben11, ben22, ben12, sigma13, sigma23}
 * 部分单元目前仅支持TL格式，壳单元还未支持UL格式
 */

class Shell : public StructuralElement
{
protected:
	std::vector<vec3d> lcoords; ///< 单元坐标系的局部坐标
	mat3d TEG; ///< 单元坐标系到整体坐标系的转换
	double theta; ///< 材料坐标系与单元坐标系的夹角
	bool MBCOUP; ///< 耦合项
	double offset; ///< 与参考平面即中面的距离偏置
public:
	Shell(int _n, Domain* _domain);

	virtual ~Shell()
	{
	}

	void GetNodeDofType(int inode, IntArray& answer) const override;
	MaterialMode getMaterialMode() override { return Shell3dMR; }
	int GetSpatialDimension() const override { return 3; }
	double getMaterialTheta() const override { return theta; }
	void computeTEM(mat3d& TEM) const override;
	/**
	 * \brief 计算壳单元刚度矩阵
	 * \param answer 壳单元刚度矩阵（目前已经支持几何非线性的部分）
	 * \param ctime 时间参数
	 */
	void ComputeStiffnessMatrix(matrix& answer, TimeStep* ctime) override;

	/**
	 * \brief 计算材料本构
	 * \param answer 计算本构方程
	 * \param rmode 材料属性用以确定那个本构
	 * \param ctime 时间参数
	 * \param ip 积分点用以获取材料点上携带的应力应变信息
	 */
	void ConstitutiveMatrix(matrix& answer, MaterialStiffnessMode rmode, TimeStep* ctime,
	                        IntegrationPoint* ip);

	void computeGeometryStiffnessMatrix(matrix& answer, TimeStep* ctime) override;

	/**
	 * \brief 小变形下用于更新应力应变
	 */
	void UpdateInternalState(TimeStep* ctime) override;

	void deflectderivate(const vec2d* G, const DoubleArray& disp, double& dwdx, double& dwdy);

	void Update(TimeStep* ctime) override;

	/**
	 * \brief 计算积分点处雅克比的值
	 */
	double det(int n, TimeStep* ctime);

	/**
	 * \brief 计算内力 （Bmat + BNL) * Dmat * (Bmat + 0.5 * BNL) * localdisp
	 * \param answer 内力
	 * \param ctime 时间变量
	 */
	void ComputeInternalForces(DoubleArray& answer, TimeStep* ctime, int userUpdate = 0) override;

	/**
	 * \brief 计算壳单元的局部坐标系(必须在子类中进行重载)
	 */
	void computeLocalSystems() override = 0;

	/**
	 * \brief 计算质量矩阵，目前对壳单元重载为集中质量矩阵
	 * \param answer 质量矩阵
	 * \param ctime 时间参数
	 */
	void ComputeMassMatrix(matrix& answer, TimeStep* ctime) override
	{
		this->ComputeLumpedMassMtrx(answer, ctime);
	}
	/**
	 * \brief 忽略旋转自由度时的质量
	 * \param answer 集中质量矩阵
	 * \param ctime 时间参数
	 */
	void ComputeLumpedMassMtrx(matrix& answer, TimeStep* ctime) override;

	void ComputeBoundarySurfaceLoadVectorAt(DoubleArray& answer, SurfaceLoad* load, int boundary, CharType type, ValueType mode, TimeStep* ctime) override;

	/*
	 * 输出积分点上的变量
	 */
	int getIPValue(DoubleArray& answer, IntegrationPoint* ip, InternalValueStateType type, TimeStep* ctime) override;

	// 根据传入的节点矢量判断是第几节点并获取其相关的转换矩阵
	int getEdgeTransformationMatrix(const IntArray& edge, DoubleMatrix& m1, DoubleMatrix& m2, IntArray& index) override;

	double* GaussWeights() const override { return &((FESurfaceElementTraits*)(m_pT))->gw[0]; } // weights of integration points
	double* Gr(int n) const { return ((FESurfaceElementTraits*)(m_pT))->Gr[n]; } // shape function derivative to r
	double* Gs(int n) const { return ((FESurfaceElementTraits*)(m_pT))->Gs[n]; } // shape function derivative to s
	double* GH(int n) const { return ((FESurfaceElementTraits*)(m_pT))->m_H[n]; }

	double r(int n) const { return ((FESurfaceElementTraits*)(m_pT))->gr[n]; } // integration point coordinate r
	double s(int n) const { return ((FESurfaceElementTraits*)(m_pT))->gs[n]; } // integration point coordinate s

	virtual void shearPart(matrix* S, int n) = 0;

	MXPOSTELEMENT getElementMxdbType() override { return SHELL_MXPOST; }

	bool NodalAveragingRecoveryModel_Extrapolation(matrix& answer, const IntArray& intvalsToExport, TimeStep* ctime, const IntArray& defaultSizes) override;
private:
	// 获取某一坐标系下材料矩阵
	// TXM:某一坐标系到材料坐标系的转换矩阵
	void computeConstitutiveAtSpecificLocalSystem(mat3d& G1, mat3d& G2, mat2d& G3,
	                                              mat3d& G4, const mat3d& TXM, TimeStep* ctime);

	// 计算形函数偏导
	double shapeGradient(int n, vec2d* G, TimeStep* ctime);

	// 计算雅克比矩阵的逆
	double invJact(int n, TimeStep* ctime, double Ji[2][2]);
};
#endif
